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Abstract 

The dynamically triangulated random surface (DTRS) approach to Euclidean quan- 
tum gravity in two dimensions is considered for the case of the elemental building 
blocks being quadrangles instead of the usually used triangles. The well-known al- 
gorithmic tools for treating dynamical triangulations in a Monte Carlo simulation 
are adapted to the problem of these dynamical quadrangulations. The thus defined 
ensemble of 4-valent graphs is appropriate for coupling to it the 6- and 8-vertex 
models of statistical mechanics. Using a series of extensive Monte Carlo simulations 
and accompanying finite-size scaling analyses, we investigate the critical behaviour 
of the 6-vertex F model coupled to the ensemble of dynamical quadrangulations 
and determine the matter related as well as the graph related critical exponents of 
the model. 
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1 Introduction 



Einstein gravity being perturbatively non-renormalizable as a field theory, con- 
structive approaches towards a quantization of gravity have been an ever more 
active field of research in the past decades 0|. The dynamical triangulations 
model in its Euclidean and Lorentzian versions has proved a successful ansatz 
for the formulation of such a consistent theory of quantum gravity 0, El- 
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Compared to the more fancy methods, such as string theory [4j and non- 
commutative geometry [5j , it is rather more minimalistic in trying to directly 
model the quantum fluctuations of space-time by a probabilistic sum over an 
ensemble of discrete, simplicial manifolds For the Euclidean case in two 
dimensions, this ensemble can be defined as the set of all gluings of equilateral 
triangles to a regular, usually closed surface of fixed topology, while counting 
each of the possible gluings with equal weight. The resulting random- surface 
model and its simplicial generalisation to higher dimensions are numerically 
tractable, for instance by Monte Carlo simulations. Furthermore, for the case 
of two dimensions the use of matrix models and generating-function techniques 
led to exact solutions for the cases of pure Euclidean gravity hi M and the 
coupling of certain kinds of matter, such as the Ising model [9, |lC| |ll| , to the 
surfaces. These two-dimensional theories generically exhibit continuous phase 
transitions on tuning the relevant coupling parameters accordingly and thus 
allow for taking the intended continuum limit. In the case of matter variables 
coupled to two-dimensional dynamical triangulations, the critical exponents 
governing the transitions are conjectured exactly from conformal field theory 
as functions of the exponents on regular lattices via the so-called KPZ/DDK 
formula [12j 



^ _ Vj - C + 24A - VT^C 
V25 - C - VT^C ' 



where A is the original scaling weight, A the scaling weight after coupling 
to gravity and C the central charge. The field-theory ansatz leading to Eq. 
(1) breaks down for central charges C > 1, an effect which has been termed 
the C = 1 "barrier", whereas the discrete model of C > 1 matter coupled 
to dynamical triangulations stays well defined. This mismatch of descriptions 
and its driving mechanism is still one of the rather poorly understood aspects 
of the dynamical triangulations model 



Ice-type or vertex models on regular lattices form one of the most general 
classes of models of statistical mechanics with discrete symmetry (for reviews 
see, e.g., Refs. 0E3)- Special cases of this class of models can be mapped 



onto more well-known problems such as Ising and Potts models or graph 
colouring problems [16]. For the case of two-dimensional lattices, several of 
these vertex models can be solved exactly, yielding a very rich and interesting 
phase diagram including various transition lines as well as critical and multi- 
critical points 0| . Thus, for two-dimensional vertex models one has the rare 
combination of a rich structure of phases and an exceptional completeness 
of the available analytical results. Hence, coupling this class of models to a 
fluctuating geometry of the dynamical triangulations type is of obvious inter- 
est, both as a prototypic model of statistical mechanics subject to annealed 
connectivity disorder and as a paradigmatic type of matter coupled to two- 
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dimensional Euclidean quantum gravity. Recently, the use of matrix model 
methods led to a solution of the thermodynamic limit of a special 6-vertex 
model, the F model, coupled to planar 4 graphs jl7;]. It was found to corre- 
spond to a C = 1 conformal field theory, i.e., it lies on the boundary to the 
region C > 1, where the KPZ/DDK solution [l2j breaks down. Also, a special 
slice of the 8-vertex model could be analysed via transformation to a matrix 
model 0]. A generalisation of this result to the gen eral parameter space of 
the 8-vertex model is currently being attempted 19, 20]. However, owed to 
the method of matrix integrals, these studies neither reveal the behaviour of 
the matter related observables and the details of the occurring phase tran- 
sitions nor the fractal properties of the graphs such as, e.g., their Hausdorff 
dimension. Especially, for the case of the 6-vertex model, which turns out to 
exhibit a phase transition of the Berezinskii-Kosterlitz-Thouless (BKT) type, 
quantities related to the staggered, anti-ferroelectric order parameter cannot 
be easily constructed, such that a detailed numerical analysis of the problem 
seems valuable. Numerically it is found here that, due to the combined effect 
of the presence of logarithmic corrections to scaling expected for a C — 1 
theory and the comparative smallness of the effective linear extent of the ac- 
cessible graph sizes, the leading scaling behaviour is obscured by extremely 
strong finite-size corrections. Thus, a very careful scaling analysis incorporat- 
ing the various correction terms has to be performed in order to disentangle 
the corrections from the asymptotic scaling form. 

Since the 6- and 8-vertex models of statistical mechanics are defined on a 
lattice with four-valent vertices, instead of considering dynamical triangula- 
tions or the dual planar, "fat" (i.e., orientable) 3 graphs, one has to use an 
ensemble of dynamical quadrangulations or the dual 4 Feynman diagrams 
as the geometry to model the coupling of vertex models to quantum gravity. 
This can be rather easily done within the framework of matrix model methods 
@ HH- For Monte Carlo studies, however, it turns out that the well established 
simulation techniques for dynamical triangulations 0, 113, are quite cum- 
bersome to adapt to the case of four-valent graphs which, therefore, only very 
scarcely have been considered in the literature 0,0- Especially, ergodicity 
for the selected set of moves has to be ensured and a method of coping with 
the observed severe critical slowing down of the dynamics, such as an adaption 
of the "baby- universe surgery" method 23[ lil , has to be devised. The details 



of these modifications to the simulation scheme will be presented in a separate 
publication |27j ]. 

The rest of this paper is organised as follows. In Sec. 2 we first review the 
basic properties of vertex models on regular lattices. We shortly discuss the 
matrix-model solution of the 6-vertex model and elaborate on the necessary 
conceptual and simulational modifications for considering vertex models on 
random graphs. Section 3 is devoted to an in-depth investigation of the BKT 
phase transition of the 6-vertex F model coupled to planar, "fat" <j) 4 graphs 
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by means of an extensive series of Monte Carlo simulations. In Sec. 4 we 
present our numerical results for the geometrical properties of the coupled 
system, such as the string susceptibility exponent and the internal Hausdorff 
dimension. Finally, Sec. 5 contains our conclusions. 



2 Vertex models on random graphs 



2.1 Vertex models on regular lattices 



An ice-type or vertex model was first proposed by Pauling 28( as a model for 
(type I) ice. In this model, the two possible positions of the hydrogen atoms 
on the bonds of the crystal formed by the oxygens, if symbolised by arrows, 
lead to six different allowed configurations around a vertex provided that the 
experimentally observed ice rule is satisfied, stating that each vertex has two 
incoming and two outgoing arrows, see, e.g., Ref. |l6j. While for the origi- 
nal ice model all vertex configurations were counted with equal probability, 
for the general 6-vertex model vertex energies e« are introduced, resulting in 
Boltzmann factors Ui = exp(— ei/ksT), where T denotes temperature and ks 
is the Boltzmann constant. Some symmetry relations are commonly assumed 
between the weights uf, in particular, given the interpretation of the arrows as 
electrical dipoles, in the absence of an external electric field the partition func- 
tion should be invariant under a simultaneous reversal of all arrows, leading to 
the identities a = u\ = uj 2 , b = 003 = U4 and c = U5 = lo§. An especially sym- 
metric version of the model assumes e a — e& = 1, e c = or a — b, c — 1. This 
so-called F model |29| turns out to exhibit an aniz-ferroelectrically ordered 
ground state. The square-lattice, zero-field 6-vertex model has been solved 
exactly in the thermodynamic limit by means of a transfer matrix technique 



(Bethe ansatz) by Lieb [30] and Sutherland j3jj. The analytic structure of the 



free energy is most conveniently parameterised in terms of the variable 0] 

A = + (2) 
lab y ' 



The free energy takes a different analytic form depending on whether A < — 1, 
— 1 <A< 1 or A > 1. Thus, phase transitions occur, whenever |A| = 1. 
The case A > 1 corresponds to two symmetry-related ferroelectrically ordered 
phases termed I and II, A < — 1 denotes an anti- ferroelectrically ordered phase 
IV and — 1 < A < 1 is attained in the disordered phase III. The latter phase 
has the peculiarity of having an infinite correlation length throughout, which 
can be traced back to the fact that it corresponds to a critical surface of the 
more general 8- vertex model 16|| . From Eq. (2) it is obvious that the F model 
exhibits a phase transition on cooling down from the infinite-temperature 



4 



point a = b = c = 1 contained in the disordered phase III to somewhere in the 
anti-ferroelectrically ordered phase IV. The transitions I — > III and II — > III 
are first-order phase transitions 0|. The transition III — > IV of the F model, 
on the other hand, exhibits an essential singularity of the free energy known 
as the BKT phase transition [32^ . 

While the ferroelectrically ordered phases exhibit an overall polarisation which 
can be used as an order parameter for the corresponding transition, the anti- 
ferroelectric order of phase IV is accompanied by a staggered polarisation 
with respect to a sub-lattice decomposition of the square lattice. That is, 
when decomposing the square lattice into two new square lattices tilted by 
7r/4 against the original one, the anti-ferroelectric ground states correspond 
to a ferroelectric ordering of the vertices of the sub-lattices with opposite 
signs of the overall polarisation of the sub-lattices. An order parameter for 
the corresponding transition can be defined by introducing overlap variables 
<Tj for each vertex of the lattice such that <7j = v $ * v®, where Vi denotes the 
arrow configuration at vertex i, v® one of the two anti-ferroelectric ground- 
state configurations and the product "*" symbolises the overlap given by 

v*v' = Y,Mv)Mv')> (3) 

k=l 



where k numbers the four edges around each vertex and Aj-{y) should be 
+ 1 or —1 depending on whether the corresponding arrow of v points out of 
the vertex or into it . Then, the spontaneous staggered polarisation P = 
(o"j)/2 = (a)/ 2 vanishes in the disordered phase and approaches unity in the 
thermodynamic limit for low temperatures in phase IV and can thus be used 
as an order parameter for the anti-ferroelectric transition. 

Vertex models on regular lattices are closely linked with different series of 
integrable models, which in turn are related to an exhaustive enumeration of 
certain conformal field theories. In fact, it turns out that the 6- vertex model, 
being the critical version of the 8- vertex model, includes in suitable generalisa- 
tions the critical points of all of the well-known two-dimensional lattice mod- 
els of statistical mechanics, including the Ising and Potts models as the most 
orominent examples. Especially, the restricted solid-on-solid (RSOS) models 



33 1, which realise each central charge of the unitary series of minimal models 
|, have been shown to asymptotically map onto the 8- vertex model, such 
that the critical RSOS models correspond to 6-vertex models. Furthermore, 
an impressive series of models in two dimensions can be mapped onto the 
Coulomb gas [i^. In these mappings, an intermediate step is always given by 
models of the solid-on-solid (SOS) type, which again can be related to vertex 
models [36|. Combining these methods, the 6-vertex model can be described 
as the common element among critical systems in two dimensions |37| . 
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2.2 Vertex models on random lattices 



Putting a vertex model onto a random four-valent graph such as the quan- 
tum gravity c/> 4 graphs imposes an additional restriction on the class of vertex 
weights that can be sensibly considered. The ferroelectrically ordered phases 
I and II of the 6-vertex model and the order parameter describing the cor- 
responding phase transition depend on the existence of a global notion of 
direction. On a random graph, this notion is maldefined. The only local orien- 
tational structure available is that of the vertices and faces of the graph. Thus, 
for an 8-vertex model coupled to quantum-gravity 4 random graphs, one has 
to assume that a = b, while the other vertex types can still be distinguished 
with only a cyclic ordering of the links around each vertex. For the 6-vertex 
model this leaves only two fundamentally different choices of models to be 
sensibly considered: the F model with e a — e& = 1, e c = and the so-called 
inverse F (IF) model with e a = ej, = —1, e c = 0, which, however, is not of 
much interest here due to its lack of an ordered phase. 

For the square lattice an order parameter for the anti-ferroelectric transition 
of the F model could be defined by a suitably calculated overlap between the 
actual state and one of the two anti-ferroelectrically ordered ground states. On 
a random graph, the corresponding ground states are not so easily found and, 
moreover, vary between different realisations of the connectivity of the graph. 
Hence, to define an anti-ferroelectric order parameter for the random graph 
case, a different and more suitable representation of the vertex model has to 
be sought. Above, the anti-ferroelectrically ordered state has been described 
as mutually opposite ferroelectric order on two complementary sub-lattices. A 
decomposition of the square lattice of this kind corresponds to a bipartition or 
two-colouring of its sites. Unfortunately, the considered random c/> 4 graphs are 
not bipartite in general, preventing an immediate application of this prescrip- 
tion. When interpreting the vertex-model arrows as a discrete vector field on 
the lattice, the ice rule for the 6-vertex model translates to a zero-divergence 
condition for this field. We thus transform the vertex model from its inter- 
pretation as a field on the links of the original lattice to a representation of 
the curl of this field on the faces of the lattice or, equivalently, the sites of 
the dual lattice. Following Stokes' theorem, this is done by integrating the 
vertex model arrows around the elementary plaquettes. By convention, pla- 
quettes are traversed counter-clockwise, adding +1 for each arrow pointing in 
the direction of motion and —1 otherwise. On the square lattice the resulting 
"spins" (or "heights") on the plaquettes can assume the values 0, ±2, ±4. This 
is demonstrated in Fig. 1. In this way, the 6-vertex model can be transformed 
to a sort of "spin model" on the dual of the original lattice. Note, however, 
that one has rather involved restrictions for the "spin" values allowed between 
neighbouring plaquettes, which would lead to quite cumbersome interaction 
terms when trying to write down a Hamiltonian. 
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Fig. 1. Transformation of the square-lattice 6- vertex model to a "spin" model 
on the dual lattice. The four links of each plaquette of the lattice are traversed 
counter-clockwise. The "spin" values written in the centres of the plaquettes are 
oriented sums of ±1 around the plaquettes. Thus, the occurring "spin" values are 
0,±2,±4. 

In the new representation, the anti-ferroelectrically ordered state of the model 
again has a sub-lattice structure. However, in contrast to the sub-lattice de- 
composition of the original representation, now the dual lattice is broken down 
into "black" and "white" sub-lattices, such that no two plaquettes of the same 
colour share a link. Then, an order parameter for the anti-ferroelectric transi- 
tion can be defined as the thermal average of the sum of the plaquette "spins" , 
e.g., on the "black" plaquettes. Reflecting the construction of the plaquette 
"spins" in Fig. 1 it is obvious that this definition of the order parameter exactly 
coincides with the original definition of Sec. 2.1 on the level of configurations. 
The difference is, however, that the new definition can be easily generalised 
to the case of arbitrary lattices, as long as their duals are bipartite. This is 
the case for the planar random 4 graphs we are considering since any planar 
quadrangulation is bipartite. Thus, we can introduce a two-colouring of the 
faces of the graphs. While for the square lattice the numbers of black and 
white plaquettes are always the same, the black and white faces of the <p A 
random graphs not necessarily occur at equal proportions. Thus, one should 
take the "spins" of both types of faces into account, however "weighted" with 
the colour of the faces. Therefore, the configurational value of the staggered 
polarisation of the F model on a planar <f) A random graph Q can be defined 
as P = | Y^vevtg*) C v S Vy where Q* denotes the dual of the graph, i.e. the 
quadrangulation, V{Q*) the set of vertices of Q*, C v = ±1 the "colour" of 
the plaquette of Q corresponding to the vertex v of Q* and S v the plaquette 
"spin" at v. Recalling the construction of the plaquette "spins", this can also 
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be written in terms of the 4 graph Q as 

p = \ E E^C/), (4) 

/GF(5) l f ef 



where F(Q) denotes the set of faces of Q, If the links of face f,Cf — ±1 
the "colour" of / and A(lf) = ±1 the direction of the vertex-model arrow 
on link If with respect to the prescribed anti-clockwise traversal of the faces. 
The thermal average (P)/2 is now taken as the order parameter of a possibly 
occurring anti-ferroelectric phase transition of the F model coupled to pla- 
nar <f) 4 random graphs. Note, however, that due to the overall arrow reversal 
symmetry of the vertex model the expectation value (P) will vanish at any 
temperature for a finite graph. Thus, for finite graphs we consider the modulus 
(|P|) instead, analogous to the usual treatment of the magnetisation of the 
Ising model. 

As mentioned above in the Introduction, a matrix model related to the F 
model coupled to planar 4 random graphs could be solved exactly in the 
thermodynamic limit ^3]. The solution is related to a transformation of the 
F model to a model of close-packed loops by using "breakups" of the vertices, 
i.e., prescriptions for connecting incoming and outgoing arrows. The original 
weights of the 6-vertex model translate into weights for the oriented loops 
by assigning a phase factor exp(z/i7r/2) to each left turn and a phase factor 
exp(— ifin/2) to each right turn of an oriented loop jTfUiij]. Here, the coupling 
fi is related to the weights of the F model as 2 

a/c = b/c= [2 cos(tt fx)]- 1 . (5) 



On the square lattice the phase factors around each loop always multiply up 
to a total of exp(±i/i27r) due to the absence of curvature. On a random graph, 
however, a loop I in general receives a non-trivial weight exp[iLiT(l)] with T(7) 
denoting the integral of the geodesic curvature along the curve I, i.e., 

71 

r(/) = — (# left turns — # right turns) . (6) 
2 



This loop expansion is related to the well-known loop representation of the 
O(n) model of Ref. [i^. There, on a regular lattice, due to the absence of 
curvature all loops receive the same constant fugacity n = 2 exp(±z/z27r), 
leading to the critical 0(n) model. On the considered random graphs this 

2 Note that, in terms of the parameter A of Eq. (2), this choice of weights covers 
only the range — 1 < A < 1, which corresponds to the disordered phase of the 
square-lattice F model. 
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picture only remains valid for the limiting case fi = 0, where the curvature 
dependence cancels. Thus, the // = point of the F model on random planar 
/graph, is equivalent to the critical 6(2) loop model 0HI |4l| and thus, by 
universality, the critical XY model 3 . Note that this corresponds to the same 
critical point a/c = b/c = 1/2 as on the regular square lattice, which is natural 
since the symmetry breaking is induced by the choice of the vertex weights. By 
means of the mentioned matrix model techniques it is found that the F model 
coupled to planar, "fat" c/> 4 graphs has a critical point for each value of the 
coupling n (corresponding to the disordered phase III), in agreement with the 
behaviour on the square lattice. Exploring the vicinity of this critical point, 
it is found that the string susceptibility exponent 7^ = for all /i, leading to 
only logarithmic divergences of the free energy ^tJ. This behaviour is indeed 
expected from the C — > 1 limit of the KPZ/DDK prediction Eq. (1). Thus, the 
general phase structure of the F model coupled to planar random 4 graphs in 
the grand-canonical ensemble of a varying number of vertices has been found 
in Ref. The existence of a BKT type phase transition at /x = was obvious 
beforehand from the equivalence to the 0(2) loop model at this point. Details 
of the behaviour of matter- related observables close to the critical point, such 
as the scaling of the staggered anti-ferroelectric polarisability, however, could 
naturally not be extracted from the matrix model ansatz. 



3 The anti-ferroelectric phase transition 



The critical point of the F model on the square lattice provided the first ex- 
ample of an infinite-order phase transition of the BKT type. By virtue of the 
loop expansion sketched above, this behaviour is expected to persist as the 
model is coupled to a random lattice. In the vicinity of a phase transition of 
this type, the usual thermal and finite-size scaling (FSS) relations are pro- 
foundly changed. Using an elaborate set of simulational techniques specially 
tailored for simulations of this model, we present a detailed scaling analysis 
of its thermal properties. As a guideline for the rather involved analysis we 
used our newly performed set of loop-cluster update simulations of the square- 
lattice F model [43], which is computationally less demanding such that much 
larger system sizes could be investigated. For a general discussion of scaling 
and FSS at an infinite-order phase transition of the BKT type we refer 
the reader to this study and references found therein 43|| . Due to the nature of 



the occurring singularities the main strengths of FSS are found not to apply to 
the BKT phase transition, and the focus of numerical analyses of the XY and 



3 Note that the loops occurring in the expansion of the O(n) model are not in 
general close packed on the lattice as are the loops of the presented loop expansion 
of the F model. However, the critical 0(2) model lies at the boundary of the dense 
phase of the 0(n) model, where loops are close packed |42T |. 
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related models has been on thermal scaling, see, e.g., Ref. |44j . In addition, 
renormalization group analyses predict logarithmic corrections to the leading 
scaling behaviour j45(, as expected for a C = 1 theory, which have been found 
exceptionally hard to reproduce numerically due to the presence of higher 
order corrections of comparable magnitude |46j . Comparing the phase transi- 
tions in the two-dimensional planar and the six-vertex F models, one should 
keep in mind that due to the dual relation of both models, the roles of high- 
and low-temperature phases are exchanged in that the F model has a critical 
low-temperature phase, whereas the high-temperature phase is massless in the 
XY model. In contrast to the XY model, the low-temperature phase of the 
F model exhibits a non-vanishing order parameter, given by the spontaneous 
polarisation of Eq. (4), such that, although the critical points of both models 
are equivalent, the magnetisation of the planar model does not correspond to 
the polarisation of the F model |43l ]. 



3.1 Simulation techniques 



For Monte Carlo simulations of two-dimensional combinatorial dynamical tri- 
angulations or the dual regular <p 3 graphs, an ergodic set of updates for simu- 
lations of a fixed number of polygons or graph vertices (canonical ensemble) is 
given by the so-called Pachner moves |47[. An adaption of the link-flip move 
for canonical simulations of triangulations to the case of quadr angulations has 
been proposed in Refs. |24],|25j]. By the construction of counter-examples it can 
be shown that the link-flip moves of Refs. 24 , 25| do not in general constitute 



an ergodic dynamics for canonical simulations of dynamical quadrangulations. 
Introducing a second type of link-flip moves, we construct an algorithm for 
canonical simulations of dynamical quadrangulations, which does not show any 
signs of ergodicity breaking (2?], @, Hi|. A scaling analysis of the thus con- 
structed dynamics reveals that its performance — as expected from a local 
algorithm — is limited by the effect of critical slowing down. To alleviate this 
problem, we adapt the non-local "baby-universe surgery" method proposed in 
Ref. [23] for triangulations to the case of quadrangulations and investigate its 



dynamical properties by means of a scaling analysis 27, 49]. For the vertex 
model part, we also employ a non-local, cluster algorithm known as "loop- 
cluster algorithm" , which is known to drastically reduce autocorrelation times 
for vertex models on the square lattice [5(j. For the application of this simu- 
lation scheme to random-lattice models, certain modifications are necessary. 
The mentioned algorithmic developments for the graph and the vertex model 
part as well as the technical details of the necessary simulational set-up will 
be discussed in a separate publication |271 |. 

For the FSS study to be presented below, we simulated a series of spherical 4 
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graphs of sizes ranging from N2 = 256 up to N2 = 65 536 vertices 4 . The sim- 
ulations where performed at several computing facilities using about 100 000 
hours of CPU time in total. 



3.2 Scaling analysis 



We assume a parameterisation of the F model coupling parameters which 
involves a temperature variable and thus sticks more closely to the language 
of statistical mechanics than to that of field theory. It hence differs from the 
parameterisation (5) used in the context of the matrix model solution, which 
only covers the critical disordered phase of the F model. Choosing the vertex 
energies as e a = e& = 1, we have a = b = e _/3 , c = 1, where (3 = l/k B T, such 
that the BKT point occurs for (3 C = In 2, both for the square-lattice model 
and, conjectured by the matrix model solution discussed in Sec. 2.2, for the 
F model coupled to planar 4 random graphs. 



3.2.1 The specific heat 

The specific heat C v of the F model coupled to planar 4 random graphs 
exhibits a broad peak of around C v ~ 0.45, shifted away from the critical 
point into the low-temperature phase 5 to a centre of (3 ~ 1.0. The peak does 
not depend on the lattice size up to very small finite-size corrections, i.e. 



15. 4 



no FSS is observed. The expected essential, non-divergent singularity 
cannot in general be resolved, since it is covered by the presence of non-singular 
background terms. This non-scaling behaviour of the specific heat is commonly 
considered as a first good indicator for the presence of an infinite-order phase 
transition 15111. 



3.2.2 Location of the critical point 

The critical coupling can be determined from the scaling of the shifts of suit- 
ably defined pseudo-critical couplings on finite graphs, see Ref. H. Here, we 



use the locations f3 x of the maxima of the staggered anti-ferroelectric polar- 
isability, defined from the generalised polarisation of Eq. (4). In terms of the 



4 The use of the variable iV~2 for this number has its origin in the general notation for 
simplicial manifolds, where denotes the number of d-simplices of the simplicial 
complex, see, e.g., Ref. 0. 

5 Note that the specific heat of the 2D XY model exhibits a peak in the high- 
temperature phase, as expected from duality. 
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inverse temperature (3 to first order one has at a BKT transition |43|, 

(3 x (N 2 )=(3 c + A p (\nN 2 )-^ } (7) 



where N 2 is the size of the graphs and p = 1/2 for the regular XY and F 
models [3 El| • For the determination of the peak positions we make use 



of the temperature- reweighting technique |52j. Note that the quoted errors 
do not cover the potential bias induced by the reweighting procedure. We 
performed simulations for graph sizes between N 2 = 256 and N 2 = 25 000 
sites, taking some 10 6 measurements after the systems had been equilibrated. 
Measurements were taken after every tenth sweep of the combined link-flip 
and "baby-universe surgery" dynamics, using "regular" graphs without self- 
energy and tadpole insertions |27|. All statistical errors were determined by a 
combined binning/ jackknife technique, cf. Ref. |53| . 

Comparing the estimated peak locations to the corresponding results for the 
square-lattice model |43| one notes that the accessible part of the scaling 
regime is strongly shifted towards lower temperatures, being rather far away 
from the conjectured critical coupling (3 C = In 2 « 0.693, cf. the "regular 
ensemble" data of Fig. 2 below. We start with fits of the simple form Eq. (7) 
without including any correction terms. Additionally, we assume p = 1/2 here 
as in the square-lattice case, which has to be justified a posteriori by the 
thermal scaling analysis. Within this scheme, the influence of correction terms 
is taken into account by successively omitting lattice sizes from the small- N 2 
side. Due to the strong corrections present, however, no fits with satisfactory 
fit quality can be found in this way such that it appears mandatory to include 
correction terms. Since the exact form of the present scaling corrections is 
not known, an effective description has to be employed. One possible ansatz 
is to relax the constraint p = 1/2, introducing p ^ p as an additional fit 
parameter. Even for this type of fit, acceptable fit qualities can only be attained 
by dropping many of the smaller graph sizes, thus strongly increasing the 
uncertainty in the estimated parameters. Additionally, we find that the fit 
results for small minimum included graph sizes iV" 2im i n partly depend on the 
choice of the starting values for the fit parameters, i.e., that the fit routine 
gets stuck in local minima of the x 2 distribution. For N 2im i a = 4096 we arrive 
at an estimate (3 C = 0.83(58), Ap = 1.7(62), and 1/p = 1.0(31) with a quality 
of Q = 0.69. Statistically, this is in agreement with the expected value j3 c = 
In 2 « 0.693 for the critical coupling, but due to the large statistical error the 
estimate is of limited significance. The result for the exponent p cannot be 
taken as a serious estimate for p, since it incorporates corrections effectively. 

For the square-lattice case, from the exact solution the leading corrections to 
the form (7) with p = 1/2 could be expressed as a power series in l/\nN 2 
0, 

P X (N 2 ) = (3 C + Ap(lnN 2 y 2 + B/ 3 QiiN 2 )~ 3 + C^ln iV 2 )- 4 , (8) 
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Table 1 

Parameter results of linear fits of the form (8) to the simulation data for the peak lo- 
cations of the staggered polarisability. Values of parameters held fixed are indicated 
by square brackets. 



^2,min 




4> 


Bp 




Q 


256 


0.8999(71) 


17.4(11) 


-69.4(47) 


[0] 


0.02 


512 


0.876(13) 


21.9(22) 


-92.0(108) 


[0] 


0.08 


1024 


0.817(24) 


33.7(46) 


-155.3(243) 


[0] 


0.72 


256 


0.779(39) 


55.4(122) 


-424.9(1140) 


918(294) 


0.32 


512 


0.693(75) 


87.0(263) 


-748.1(2647) 


1838(741) 


0.39 



hence we consider this form for the random graph data here as well. As can 
be seen from the collection of fit parameters in Table 1, this form provides 
a good description of the data, although some of the statistical errors of the 
fit parameters become very large. Neglecting the second correction first, i.e., 
holding C/3 = fixed, the results are stable on successively omitting data 
points from the small- N 2 side, and the resulting estimates for the transition 
temperature are slowly drifting towards the asymptotic value f3 c = In 2 = 
0.693. . .. Nevertheless, the result for, e.g., N 2 , min = 1024, (3 C = 0.817(24) is 
still far from being compatible with the asymptotic result in terms of the 
statistical error. Including the fourth-order term of (8), on the other hand, 
further reduces the estimates for f3 c to the extent of being compatible with 
f3 c = In 2, however at the price of largely increased statistical errors. For 
-^2,min > 512, the fits get very unstable, such that we quote as our final result 
from this approach f3 c = 0.693(75) for A^ 2 , m in = 512. If we finally fix f3 c at 
its asymptotic value, for Cp = we reach a fit quality of Q = 0.01 only 
at A^min = 2048, while with variable Cp, Q = 0.52 is reached already at 
-^2,min = 512. This clearly shows that both correction terms are necessary for 
resolving the scaling corrections, but the accuracy of the present data is only 
marginally sufficient to do so. It should be noted that also the other types of 
fits presented here still yield good quality-of-fits when fixing the parameter f3 c 
at In 2. For example, a fit of the form (7) with variable exponent p to the data 
with iV 2 , m in = 2048 gives Ap = 1.071(81), 1/p = 0.541(35), and Q = 0.84. 



3.2.3 Universality of the critical coupling 

One might be tempted to suspect that the observed rather large deviations of 
the finite-size positions of the polarisability maxima from the expected value 
(3 C = In 2 « 0.693 are due to the fact that we use graphs of the regular ensem- 
ble, i.e. those without self-energy and tadpole insertions, whereas the matrix 



model calculations of Ref. |17| naturally concern graphs of the unrestricted 
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singular ensemble. Indeed, quite generally one does not expect the critical 
coupling of a model to be universal. In particular, for the Ising model coupled 
to dynamical polygonifications or the dual graphs, the location of the observed 
transition does depend on whether one considers spins located on the vertices 
of triangulations, quadrangulations, 3 or <p A graphs 24 , 25| . Additionally, 



depending on the considered ensemble of graphs with respect to the inclusion 
or exclusion of certain types of sin gula r contributions, one arrives at different 
values for the critical coupling 1(J 11 , 54, HI]. However, the situation is quite 



different for the case of the F model coupled to random lattices. As has been 
mentioned above in Sec. 2.2, in the matrix model description of the problem, 
the matrix potential becomes equivalent to that of the 0(2) model in the limit 



jj, = |17| . which corresponds to the choice a/c = b/c = 1/2 or (3 C = In 2. 
Thus, renormalizing the matrix model to remove some or all of the singular 
graph contributions does not change the location of the BKT point. 

We have not performed extensive simulations of graphs of the "singular" en- 
semble including self-energy and tadpole insertions to demonstrate this be- 
haviour numerically. This is due to the fact that simulations for graphs of the 
singular ensemble are by orders of magnitude less efficient for the considered 
graph sizes than simulations of the other graph ensembles due to details of the 
implementation of the simulation scheme, cf. Ref. |22j]. Nevertheless, we car- 
ried out some simulations for smaller graph sizes and analysed the FSS of the 
peak locations of the staggered polarisability just as for the case of "regular" 
graphs. The corresponding FSS data are shown in Fig. 2 together with the re- 
sults for regular graphs. Using Eq. (8) with Cp = 0, a fit to the data including 
all five points from N 2 = 128 to N 2 = 2048 yields the estimate p c = 1.01(11), 
Q = 0.92, letting Cp vary gives (3 C = 0.83(69), Q = 0.76, which is in principle 
in agreement with (3 C = In 2, although very inaccurate. Note that from Fig. 
2 the finite-size corrections for the singular graph case are much larger than 
those for the regular graph model. This is in contrast to previous observations 
for the case of the Potts model coupled to random triangulations |56| and 
the resulting common belief that the inclusion of singular graph contributions 
generically reduces FSS corrections. 

As has been mentioned in the Introduction, the reason for the observed very 
slow approach to the expected asymptotic behaviour lies in the double effect 
of the presence of logarithmic corrections to scaling and the small effective 
linear extent of the highly fractal lattices. In principle it should be possible to 
resolve the resulting scaling corrections by including higher-order correction 
terms in the fit ansatze. However, it must be admitted that, refraining from 
any artificial "good-will" tinkering with the fit parameters, the accuracy of the 
present data is not sufficient for reliable many-parameter, possibly non-linear 
fits. The strength of this combined effect is nicely demonstrated numerically 
by the fact that the fits to the FSS of the polarisability peak locations with f3 c 
fixed to its true value (3 C = In 2 come as close as /3 X (N 2 ) = 0.7 to the critical 
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Fig. 2. Finite-size approach of the peak locations of the staggered polarisability of 
the F model on c/> 4 random graphs with ("singular ensemble") and without ("reg- 
ular ensemble") tadpole and self-energy insertions. The solid lines show fits of the 
functional form (8) to the data. 

value only for graph sizes N 2 ~ 10 100 for the form (8) with variable Cp or 
even N 2 ~ 10 5000 for the form (7) with variable exponent p. Instead of figuring 
out more elaborate fits, we try to disentangle the two correction effects by a 
comparison to the square-lattice model, where only the logarithmic corrections 
are present, but the considered lattices are not fractal |43[. For this purpose, 
we plot in Fig. 3 the polarisability peak locations as a function of the root 
mean square extent of the considered lattices defined as 

K ' N2 \ E r ZoGn(r) 




which is the relevant measure for the linear extent of the graphs. Here, we 
take the geometrical two-point function Gu(r) as the number of graph vertices 
with a geodesic link distance r from a randomly chosen reference point po- The 
root mean square extents (r 2 ) 1//2 are related to the number of graph vertices 
according to (r 2 ) ~ N^ dh , which defines the internal Hausdorff dimension dh- 
Due to the fractal structure of the random graphs, largely differing values of 
the root mean square extent (r 2 ) 1 / 2 are found for them in comparison to square 
lattices with the same number of vertices N 2 . For the latter data (which are 
basically exact) this scaling ansatz without inclusion of any correction terms 
yields dh = 2.000(20), where the error reflects discretisation effects for small 
lattices. For the case of 4 random graphs the fit yields dh = 3.336(11). Note, 
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Fig. 3. Collapse of the FSS approach of the scaling of the peak locations of the 
staggered anti-ferroelectric polarisability of the F model on random c/> 4 graphs (left 
scale) and on the square lattice (right scale). The data for the square-lattice model 
are taken from a set of simulations presented in Ref. [4^. 

however, that the result for dh is slowly increasing as more and more of the 
small- N 2 lattices are excluded and we expect the true value of the Hausdorff 
dimension to be somewhat larger, see Refs. jH, IH, and Sec. 4 below. 
Hence, in order to obtain results for the F model at comparable linear extents 
of the square and random lattices, one has to consider rather small volumes 
for the square-lattice case. For the comparison we use L x L square lattices, 
where the edge lengths L were chosen such that the resulting root mean square 
extent comes as close as possible to the (r 2 ) 1 / 2 values for the corresponding </> 4 
random graphs. The volumes of the 4 random graphs were chosen between 
N 2 = 256 and N 2 = 8192, increasing in powers of two. 

In Fig. 3 we present a comparison of the FSS approach of the peak locations 
of the polarisability for the 4 graph and square-lattice |43| models plotted 
as a function of the linear extent (r 2 ) 1 / 2 of the lattices. Here, the abscissae of 
the plot have been scaled such as to account for the difference in the overall 
correction amplitude, but assuming the same value In 2 for the offset. From 



the two simulation points near 
amplitudes as 6 



2\l/2 



(r 2 ) 



10 we find the ratio of the correction 



A 



/3*(N 2 = 1024) - In 2 
P$(N 2 = 324) -In 2 



4.23, 



(10) 



6 These two simulation points have been chosen since there the difference in (r 2 ) 1 / 2 
between the square and random lattices is minimal within the set of considered 
lattice sizes. 
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where denotes the peak position for the random 4 graph model and 
the value for the square lattice. The thus achieved collapse of the FSS data is 
obvious from Fig. 3. Consequently, we come to the clear conclusion that the 
larger deviations of the peak locations for random graphs are simply due to an 
about four times larger overall amplitude of the correction terms as compared 
to the square-lattice model, the details of the FSS approach being otherwise 
surprisingly similar between the two considered lattice types. Especially, the 
fact that for the 4 graph case the asymptotic value j3 c = In 2 cannot be 
clearly resolved by the considered fits to the data is an obvious consequence 
of the comparative smallness of the accessible lattice sizes in terms of their 
effective linear extents (r 2 ) 1 / 2 . To underline this finding, we performed fits 
of the simple form (7) to the data for both types of lattices (there are not 
enough data points for fits with correction terms), including sizes starting 
from the points near (r 2 ) 1 / 2 « 10, which result in estimates (3 C = 0.7554(18) 
for the square lattice and (3 C = 0.9416(89) for the random graphs. In terms 
of the quoted statistical errors these are obviously both far away from the 
asymptotic result. The deviation from f3 c = In 2 is, however, just about four 
times larger for the random graph case than for the square-lattice model, in 
agreement with the previous discussion of the scaling collapse of Fig. 3. 

3.2.4 Critical energy and specific heat 

As an aside, we note that for the largest </> 4 random graphs we have simulated, 
i.e., for N2 = 65 536, at (3 = f3 c = In 2 we find the following values of the 
internal energy and specific heat per site, 

U(p = In 2) = 0.333355(11), C v {/3 = In 2) = 0.2137(12). (11) 

Comparin g th ese results to the values found analytically for the square-lattice 
F model H, U(J3 e ) = 1/3, C V (J3 C ) = 28(ln2) 2 /45 « 0.2989, we see that 
U(/3 = In 2) is very close to the value found for the square lattice, whereas 
C v {j3 = In 2) is far away from the square-lattice result. On the basis of these 
findings, we conjecture that the critical value of the internal energy of the 
F model is not affected by the coupling to random graphs, while the critical 
specific heat is. Thus, as one would expect, the critical distribution of vertex 
energies naturally changes its shape on moving from the square-lattice to the 
random graph model, but, curiously, its mean is not shifted by this procedure. 
Interestingly, this situation seems to be specific to the critical point (3 C = 
In 2 common to both models, where the two curves cross. For other inverse 
temperatures the square-lattice and random graph energies diverge, see Fig. 4. 
This probably indicates the presence of an additional symmetry at criticality. 
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Fig. 4. Temperature dependence of the internal energy U of the square-lattice and 
random <^> 4 graph F models. Simulations have been performed for a N2 = 46 2 = 2116 
square lattice and random graphs with N2 = 2048 sites. 

3.2.5 FSS of the polaris ability 

On coupling the vertex model to quantum gravity we expect a renormaliza- 
tion of the critical exponents as prescribed by the KPZ/DDK formula (1). 
In Ref. [12] KPZ/DDK focus on conformal minimal models with C < 1 cou- 
pled to the Liouville field, but their work should also marginally apply to the 
limiting case C — 1 of the model considered here. The KPZ/DDK formula 
prescribes a dressing of the conformal weights on coupling a matter system 
to the fluctuating background. To find the usual critical exponents from the 
weights, one assumes that the well-known scaling relations stay valid and thus 
arrives at 



2A f 



a: 



1-A £ 



P 



A, 



7 



2A, 



1- A, 



1-A e ' ' 1-A 6 
2 - 7] = (1 - 2A P )d h . 



(12) 



Here, A e denotes the weight of the energy operator and Ap symbolises the 
weight of the scaling operator corresponding to the spontaneous staggered 
polarisation Pq, which here takes on the role of the magnetisation operator o~ 
of magnetic models. For the special case of the infinite-order phase transition 
considered here, the usual exponents written above are not well-defined in the 
sense of describing power-law singularities. However, the corresponding FSS 
exponents, i.e., (3/dhV = Ap and ^jd^v = 1 — 2Ap, still have a well-defined 
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meaning. From the exponent 13/dhV = 1/4 for the square-lattice F model 
(with dh = d = 2) jiij], we find Ap = 1/4. Note that this weight is different 
from the weight A a = 1/16 found for the magnetisation of the critical XY 
model in two dimensions, see e.g. Ref. |6(|. For central charge C = 1 from 
Eq. (1) one arrives at Ap = 1/2 and the dressed critical exponents become 
fi/dhi 1 = Ap = 1/2 and 7/rf^z/ = 1 — 2Ap = 0, implying a merely logarithmic 
singularity of the staggered polarisability for dynamical graphs. 

For a numerical check of these conjectured exponents, there are the two prin- 
cipal possibilities of considering the FSS of the staggered polarisability at its 
maxima for the finite graphs or at the fixed asymptotic transition coupling 
f3 c = In 2. While in the asymptotic regime both approaches are expected to 
lead to identical results, this is not at all obvious in the presence of large, not 
completely controlled correction effects for the accessible graph sizes. In both 
cases, by analogy to the situation on the square lattice |43j we start from an 
FSS form including a leading effective correction term, namely, 

X (N 2 ) = A x N^(\nN 2 r-, (13) 



where xC^) is taken to be either the peak value as a function of (3 or the value 
at (3 — /3 C — In 2. We consider the peak value case first, taking the simulation 
results for the graph sizes N 2 = 256, . . . , 25 000. Omitting the correction term, 
i.e., forcing uj x = 0, and trying to control the effect of corrections to scaling by 
successively omitting data points from the small- N 2 side, results in quite poor 
fits with an exponent estimate 7 fd^v ~ 0.7 steadily decreasing with increasing 
lower cut-off A^ 2 , m in- Allowing the effective correction exponent uj x to vary, the 
resulting leading exponent estimate •y/d h i / is considerably reduced, still show- 
ing a tendency to decline as iV^min is increased, cf. Table 2(a). However, the 
fit quality is still not very good and the resulting exponent estimate for, e.g., 
-/V"2,min = 2048, ^ j d^v = 0.301(79) is not consistent in terms of the statistical 
error with the purely logarithmic singularity expected from the KPZ/DDK 
prediction. These results in principle might be improved by includin g co rrec- 
tions of the form l/(ln N 2 ) n , n = 1, 2, ... as in the square-lattice case 43], but 
the present data are not precise enough to reliably fit these terms. 

For the data at fixed coupling (3 C = In 2, simulations up to slightly larger 
graph sizes could be performed since no reweighting analysis is necessary there. 
Hence, results are available for graph sizes between N 2 = 256 and N 2 = 32 768 
sites, increasing by powers of two. For the constrained fits of the functional 
form (13) with u x — we do not find a quality-of-fit of at least 10" 2 for N 2<min 
up to 4096 and thus do not consider this form further. The parameters of fits 
including the logarithmic term are collected in Table 2(b), revealing that the 
functional form including a logarithmic correction fits the data rather well 
already for quite small values of A^min, leading to exponent estimates "f/dhV 
compatible with the conjecture 7/c^z/ = in terms of the quoted statistical 
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Table 2 

Results of fits of the functional form (13) to the simulation data for the staggered 
polarisability. (a) Fits to the data at the polarisability peak locations, (b) Fits to 
the data at the asymptotic critical coupling (3 = (3 C = In 2. 



-^2,min 


Ax 


l/dhV 


u x 


Q 


256 


0.1975(97) 


0.4749(81) 


1.698(55) 


0.00 


512 


0.116(14) 


0.406(16) 


2.22(12) 


0.00 


1024 


0.039(12) 


0.281(37) 


3.24(30) 


0.24 


2048 


0.047(37) 


0.301(79) 


3.07(68) 


0.16 



-^2,min 


Ax 






Q 


256 


0.491(19) 


0.0194(55) 


2.117(40) 


0.66 


512 


0.543(42) 


0.0304(91) 


2.026(72) 


0.91 


1024 


0.569(75) 


0.035(14) 


1.98(12) 


0.85 



errors. In fact, if we assume a purely logarithmic increase of xi^), be., if we 
fix j/dhV = 0, the data yield good-quality fits for A^min ^ 512; for iV^min = 
2048 the parameters of this purely logarithmic fit are A x = 0.3960(96), uo x = 
2.295(11), with Q = 0.39. 

The simulation data at /3 = In 2 together with this last fit are shown in Fig. 5. 
Note that for the peak-height data discussed before, such a purely logarithmic 
fit is not possible with acceptable values of Q. To enable a somewhat better 
judgement of the observed discrepancy between the scaling at the peak max- 
ima and at f3 c = In 2, we considered the same two lines for the square-lattice 
model |43|], using a range of lattice sizes comparable to that of the random 
graph case in terms of the effective linear extents as it has been discussed 
in Sec. 3.2.3. Fitting the functional form (13) with variable u x to these two 
square-lattice data sets, we find ^jdhV = 0.475(46) for the scaling at (3 = In 2 
also considered above, but an estimate of 7/0^ = 0.598(36) from the scaling 
of the peak values of x- Thus, also for the square-lattice model, the scaling 
of the peak values yields an exponent estimate lying off the expected result 
{pijdhU = 1/2 in this case), while fits at the critical coupling are in good agree- 
ment with the expectations. This is in agreement with the general observation 
of enhanced correction amplitudes of the random graph model compared to 
the square-lattice case reported in Sec. 3.2.3. 
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Fig. 5. Finite-size simulation data of the polarisability of the F model on random 
cj) 4 graphs at the asymptotic critical coupling C = In 2. The solid curve shows a fit 
of the form (13) to the data, where j/cl^u = was kept fixed. 

3.2.6 FSS of the spontaneous polarisation 

For the scaling of the spontaneous polarisation the situation is found to be 
quite similar to the above discussed case of the polarisability. We assume the 
same leading FSS form as in the square-lattice case ji^, i.e., 



where, again, Po{N-i) is taken to be either the value at the peak position of 
the polarisability or, alternatively, the result at the asymptotic critical cou- 
pling (3 C = In 2. Fits without the logarithmic correction term (up = 0) show 
unacceptable quality throughout the whole region of choices of the cut-off 
^V2,min and for both FSS series. For the polarisation at the peak locations of 
the polarisability, even fits including the logarithmic correction term of Eq. 
(14) show very poor fit quality and estimates for P/dhf which are clearly too 
small compared to the KPZ/DDK prediction (3/dhV = 1/2 in terms of their 
statistical errors. We attribute this to the generally more pronounced correc- 
tions for the values at the polarisability peak locations already noted above. 
In addition, however, the non-divergent behaviour of the polarisation makes 
it even harder to resolve the correction terms properly, and the possible pres- 
ence of systematic reweighting errors (bias) has much more severe effects here 
due to the higher statistical accuracy of the polarisation estimate. Again, the 
analogous analysis of the FSS of the square-lattice model reveals a similar 




(14) 
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Table 3 

Parameters resulting from fits of the form (14) to the finite-graph spontaneous 
polarisation at the infinite- volume critical coupling (3 C = In 2. 





M 


f3/d h v 


up 


Q 


256 


1.583(35) 


0.4633(30) 


0.726(22) 


0.74 


512 


1.658(68) 


0.4581(50) 


0.684(39) 


0.91 


1024 


1.58(11) 


0.4633(79) 


0.728(64) 


0.98 


2048 


1.48(23) 


0.469(15) 


0.779(134) 


1.00 



behaviour for comparable graph sizes in terms of the linear extent, however 
with the size of the deviations from the expected result being much smaller. 

Table 3 shows the parameters resulting from least-squares fits of Eq. (14) to 
the simulation data at the fixed coupling f3 — /3 C = ln2. The overall quality of 
the fits is much better than for the data at the polarisability peak locations 
discussed before. This is at least partially due to the fact that for the results at 
fixed coupling no bias effects induced by a reweighting procedure are present. 
We do not observe a clear overall drift of the exponent estimate (3 jdhV resulting 
from the fits as a function of the cut-off A^ 2 min and the quality-of-fit is found 
to be exceptionally high already for small values of A^ 2 , m in- The result for 
-^2,mm = 2048 is consistent with the KPZ/DDK conjecture [3/dhV = 1/2 
within about two times the quoted standard deviation. We note that the 
estimated correction exponents u x and u PQ are found to be clearly different 
from each other. In fact, from the exact solution of the square-lattice model, 
both exponents are found to be different even asymptotically j43|. In addition, 
both exponents effectively capture the presence of sub-leading corrections for 
the two observables, leading to the occurrence of further differences. 

3.2.7 Thermal scaling 

In order to extract information about the critical exponent p and possibly to 
find additional evidence for the location of the critical point, we tried to per- 
form a thermal scaling analysis and considered the dependence of the staggered 
anti-ferroelectric polarisability on the inverse temperature j3 in the vicinity of 
the critical point. Since the high-temperature phase of the F model coupled to 
4 random graphs is expected to be critical as for the case of the square-lattice 
F model, such a scaling analysis has to be performed on the low-temperature 
side of the polarisability peak. As for the square-lattice model [431 ]. we find 
scaling throughout the high-temperature phase. Due to an exponential slowing 
down of the link-flip and "baby-universe surgery" dynamics of the </> 4 graphs 
above (3 C |22j], simulations cannot proceed arbitrarily deep into the ordered 
phase. Up to the attainable inverse temperatures of about (3 = 1.4, we still 
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Fig. 6. Thermal scaling of the polar isability of the random graph F model for graphs 
with N2 = 30 000 sites. The curve shows a fit of the function (15) to the data, where 
(3 C = In 2 and p = 1/2 have been kept fixed. 

observe strong finite-size effects and no asymptotic collapse of the curves for 
different graph sizes, which is again attributed to the large fractal dimension 
of the graphs. 

The requirements of a proper thermal scaling analysis of the polarisability 
resulting from these observations are almost impossible to fulfil: one has to 
keep enough distance from the critical point for the linear extent of the graph 
to be large compared to the correlation length of the matter part to keep finite- 
size effects under control and, on the other hand, one should not proceed too 
deep into the ordered phase such as not to leave the thermal scaling region in 
the vicinity of the critical point. Thus, one would have to go to huge graph sizes 
to get rid of these constraints to a practically acceptable extent. Nevertheless, 
we attempt a thermal scaling analysis of the polarisability from simulations of 
graphs of size N2 = 30 000 with inverse temperatures ranging from (3 = 0.9 up 
to (3 — 1.6 taking about 800 000 measurements at each (3. From the square- 
lattice results one expects the scaling form jl^ . 

!**(/?) + -&)-', (15) 

which should hold for f3 — > (3^ as N 2 — > 00 and where logarithmic corrections 
have already been omitted. We find it impossible to reliably fit all four of the 
parameters involved in Eq. (15) to the available data. Varying the starting 
values we find a multitude of local minima of the \ 2 distribution, such that 
virtually any result can be "found" for (3 C and p in this way. Fixing one or the 
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other of the two parameters at the expected values (3 C = In 2 or p = 1/2, the 
fits become more stable. The dependency on the range of included values of 
f3 is found to be rather small and for (3 > 1.25 we arrive at the fit parameters 
A x = -101(4662), B x = 106(4662), p = 0.02(103), and Q = 0.03, for (3 C fixed 
at In 2 or at the parameters A x = -86(1083), B x = 324(5744), f3 c = -11(147), 
and Q = 0.04, with p fixed at 1/2. Obviously both fits are not very useful, 
such that we are finally forced to fix both parameters, f3 c and p, at their 
expected values to find A x = 0.91(41), B x = 4.20(33), and Q = 0.03. This fit 
is shown in Fig. 6 together with the simulation data. Thus, the best we can 
conclude about the thermal scaling behaviour of the polarisability is that there 
is no obvious contradiction with the expectations concerning the parameters 
(3 C and p. However, in view of the fact that already for the regular lattice 
model thermal scaling fits were not at all easily possible 0], this finding is 
probably not too astonishing. 



4 Geometrical properties 

The annealed nature of disorder applied to the vertex model via its placement 
onto dynamical 4 random graphs induces a back-reaction of the matter vari- 
ables onto the underlying geometry and thus a possible change in the (local 
and global) geometrical properties of the graphs. Since the general mechanism 
of matter back-reaction onto the graphs is the tendency to minimise interfaces 
between pure-phase regions of the matter variables, a strong coupling between 
matter and graph variables is generically only expected if the combined system 
of spin model and underlying geometry is critical. Thus, the universal graph 
properties such as the graph-related critical exponents should remain at the 
values of pure Euclidean quantum gravity, unless the coupled matter system 
has a divergent correlation length |6l| . As indicators for changes of the geom- 
etry of the coupled system, we consider the co-ordination number distribution 
as a typical local property, as well as the string susceptibility exponent and 
the Hausdorff dimension as global geometrical features. 

4-1 The co-ordination number distribution 

The distribution of co-ordination numbers of the quadrangulations, which has 
been extensively considered for the case of pure 4 graphs [49] , could be possi- 
bly altered by the back-reaction of a coupled matter model. In particular, for 
the case of the vertex model considered here, the ice rule forbids certain link- 
flip update moves and thus changes the distribution Pn 2 (q) of co-ordination 
numbers. The vertex configurations forbidden by the ice rule effectively carry 
infinite energy, such that they stay excluded even in the infinite-temperature 
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Fig. 7. Fraction re 2 of faces of length two of planar 4 "regular" random graphs with 
a coupled F model as a function of the inverse temperature ft. The drawn error bars 
are mostly covered by the size of the symbols. The solid line shows the value of n 2 
for the case of pure (ft random graphs of iV 2 = 2048 sites. 

limit ft — *■ 0. Thus, in contrast to, e.g., an Ising model a full decoupling of 
graph and matter variables for high temperatures does not occur here due to 
the entropic instead of energetic nature of the matter-graph back reaction. 

From our numerical simulations we find that on the scale of the whole distri- 
bution P/v 2 (?) no changes as a function of the inverse simulation temperature 
ft can be distinguished and the distribution looks identical to that of pure 
planar 4 graphs 0|. However, -P/v 2 (g) can be determined to high precision, 
and concentrating on a single point of the distribution, e.g., q = 2, a clear 
variation with the inverse temperature ft can be resolved, cf. Fig. 7. Also, in 
terms of the quoted statistical errors, which are of the order of 10~ 5 for the 
measurements of n 2 = P/v 2 (2), the pure graph result of n 2 = 0.296 365(32) 
for iV 2 = 2048 49] is very far away from the whole of the shown variation 
of the F model case. We find a peak of n 2 around ft ~ 0.7 with only rather 
small variations with the size of the considered graph. A similar peak of the 
fraction P/v 2 (3) of i/iree-faces for different spin models coupled to dynamical 
triangulations has been observed before, see Ref. 

Since a pronounced back-reaction of the matter variables onto the underlying 
graphs is only expected at criticality, we interpret the location of the observed 
peak of n 2 (/3) as a pseudo-critical point ft n2 which should scale to the asymp- 
totic critical coupling ft c = In 2. As for the thermal scaling analysis of Sec. 3.2, 
the precise location of the maxima can be determined from the simulation 
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data via reweighting. This has been done for the data from simulations of 
graphs of sizes between N 2 = 256 and N 2 = 4096 sites with time series of 
lengths between 8 x 10 5 and 4 x 10 6 measurements. We find only very small 
changes of this peak position on variation of the size of the graphs, such that 
within the present statistical errors j3 n2 can be considered constant. Thus, we 
do not perform a finite-size scaling fit to the data of the peak locations, but 
instead quote the result from the largest considered lattice as an estimate for 
the asymptotic critical coupling, namely 

(3 n2 = 0.6894(54), (16) 



resulting from the simulations for N 2 = 4096. This is in nice agreement with 
the expected value of j3 c = In 2 « 0.693 and almost two orders of magnitude 
more precise than the results found above from the scaling of the polarisability 
peak locations. 



4-2 The string susceptibility exponent 



In the grand-canonical ensemble of the dynamical polygonifications model 
the string susceptibility exponent 7 S governs the leading singularity of the 
partition function for spherical graphs via Z{p) ~ [p, — /i ) 2_7s 0, where \x 
denotes the chemical potential accounting for the cost of the insertion of a new 
vertex. Thus, a direct measurement of 7 S requires computationally demanding 
simulations with a varying number of polygons or graph vertices. Additionally, 
since a shift of 7 S due to the presence of some matter variables can only be 
expected at criticality, a numerical setup for the detection of such a change 
needs to tune two coupling constants, namely /i and /3, to criticality. Due to 
the combination of these two problems a reliable estimation of 7 S from grand- 



canonical Monte Carlo simulations has proved difficult, see e.g. Ref. 63 



It could be shown, however, that the string susceptibility exponent is related 



to the "baby-universe" structure of the dynamical polygonifications [64J. This 
observation can be turned into a method for the determination of 7 S from sim- 
ulations at a fixed number of polygons or graph vertices (canonical ensemble) 
[6~l| . The basic building blocks of this '"baby-universe" structure are taken as 
so-called "minimal- neck baby universes" (minBUs), which we define as sub- 
graphs which typically contain a "macroscopic" number of vertices, but are 
connected to the main graph body by only four links for the case of dynamical 
quadr angulations. A simple decomposition argument of the graphs into "baby 
universes" yields the following scaling relation for the distribution {un 2 {B)) 
of volumes B contained in minBUs of the ensemble of pure graphs of size N 2 
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0, 

(n N2 (B))~N^[B(N 2 -B)} 



(17) 



where B ^> 1 and N 2 — B 3> 1 is assumed. Also, it can be shown that the 
same relation should hold for the case of C < 1 conformal matter coupled to 
the polygonifications or dual graphs with j s then denoting the corresponding 
dressed string susceptibility exponent j§4[. For the limiting case C — 1, on the 
other hand, it is argued in Ref. (64| that the distribution of minBUs should 
acquire logarithmic corrections and look like, 

(n N2 (B)) ~ N^'[B(N 2 - B)Y"~ 2 \hyB ln(iV 2 - E)f, (18) 



with k = —2. An estimate ^jv 2 (B) for the volume distribution of minBUs can 
be easily found numerically from a decomposition of the graphs into "baby 



universes". When the minBU surgery algorithm (cf. Refs. [27j,|49() is applied, 
such an estimate can even be produced as a simple by-product of the updating 
scheme. Then, an estimate for 7 S can be found from a fit of the conjectured 
functional form (17) or (18) to the estimated distribution njy 2 (B) |6l|- In or- 
der to honour the constraints B 3> 1 and N 2 — B 3> 1 of Eqs. (17) and 
(18) one has to introduce cut-offs B min and B max , such that only data with 
-B m i n < B < B max are included in the fit. Here, the choice of the lower cut-off 
-B m in is found to be much more important for the outcome of the fit than 
the choice of B m8/X . We use the following recipe for the determination of the 
cut-offs: as a rule of thumb, we choose -B max = N 2 /8, which has turned out 
to be a good initial guess for most situations. With -B max fixed, the lower 
cut-off -Bmin is steadily increased from B min « 0, monitoring the effect of 
those increases on the resulting fit parameters, especially the estimated string 
susceptibility exponent j s . Finally, with the resulting value of -B mm fixed, a 
second adaption of B max is attempted, usually changing 5 max by factors of 
two or one half. Additionally, the quality-of-fit parameter Q is utilised as an 
indicator of whether neglected corrections to scaling are important for the 
considered window of minBU volumes B. As far as corrections to the lead- 
ing scaling behaviour are concerned, it is speculated in Ref. [6l| that a good 
effective description of the leading correction term results from the substitu- 
tion B ls ~ 2 — ► B ls ~ 2 [l + D ls /B\. Hence, the actual fits were performed to the 
functional form 

hin N2 (B) = A la + ( 7s - 2) In [B(N 2 - B)] + (19) 
for C < 1, resp. to the form 

\nn N2 {B) = A la + ( 7s - 2) \n[B(N 2 - B)] + k ln[ln B ln(N 2 - B)] + — ^, (20) 

B 
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Table 4 

Parameters of fits of (19) to the simulation data for the distribution n^ 2 (B) of 
minBUs for pure ^ random graphs. The small values of the quality-of-fit parameter 
Q for the two largest graph sizes are a side effect of the cross-correlations in n^ 2 {B). 



N 2 


Bmm 


^max 




7s 




Q 


1024 


60 


128 


18.36(49) 


-0.474(40) 


-2.9(30) 


0.79 


2048 


70 


256 


20.34(14) 


-0.495(10) 


-3.8(12) 


0.56 


4096 


70 


512 


22.030(90) 


-0.4915(63) 


-3.78(74) 


0.05 


8192 


100 


1024 


23.853(72) 


-0.4977(47) 


-4.80(87) 


0.04 



for the limiting case of C = 1. Here, the dependency on the total volume N 2 
has been condensed into the constant A la . Note that both of these fits are 
linear and the number of data points is of the order of 10 3 for the lattice sizes 
we have considered, such that a fit with four independent parameters is not 
unrealistic. In Eq. (20) we keep free parameter since the value k = — 2 

is only a conjecture and, additionally, further corrections to scaling can be 
covered in an effective way by letting k vary. 

4-2.1 Results for pure <f) A graphs 

Matrix model calculations for purejplanar dynamical triangulations yield the 
exact result 7 S = —1/2, cf. Ref. [2|. As a gauge for the method and as a 
check for the expected universality of 7 S with respect to the change from 
triangulations to quadrangulations, we apply the described technique first to 
the case of pure 4 random graphs. We adapt the lower and upper cut-offs B min 
and -B max iteratively as described above, taking into account that the usual 
error estimates of least-squares fits of (19) to the data could be misleading 
due to the apparent correlations of the points of n N2 (B) for different sizes 
B of the minBUs, which genericaliy lead to an underestimation of variances. 
We refrain from an additional extrapolation of the resulting estimates of 7 S 
towards B min — > oo suggested by the authors of Ref. [6l| since we do not see 
a proper justification for a specific extrapolation ansatz and in general find 
extrapolations of noisy (and here also strongly correlated) data questionable. 

Statistically reliable error estimates for 7 S are found by jackknifing over the 
whole fitting procedure: first the upper and lower cut-offs in B are determined 
as described using the full estimate Tin 2 (B). Then, of the order of ten jackknife 
blocks are built from the time series the estimate njsr 2 (B) is based on and fits 
with the same constant cut-offs are performed for each block to yield jackknife- 
block estimates of 7 S and the other fit parameters. Table 4 summarises the 
final results for pure </> 4 graphs of sizes N 2 = 1024 up to N 2 = 8192, taking 
about 10 9 x N 2 minBUs into account for each graph size. Obviously, finite- 
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Fig. 8. Estimates of the string susceptibility exponent 7 S from fits of (19) to the 
measured distribution of minBUs for graphs of size N 2 = 2048 coupled to the F 
model. The displayed error bars do not represent the full statistical error. The true 
amount of statistical fluctuations is indicated by the three data points with larger 
error bars at couplings (3 = 0.2, (3 = In 2, and = 1.4, where the errors have been 
evaluated by a full jackknife analysis. Note that the displayed exponent estimates in 
the high-temperature phase are effective exponents since there are large finite-size 
corrections (see text). 

size effects are relatively weak here, and we quote as final result the value for 
N 2 = 8192, 7 S = -0.4977(47), which is perfectly compatible with 7, = -1/2. 

4-2.2 Results for the F model case 

For the case of the F model coupled to the 4 graphs, we expect a variation of 
the string susceptibility exponent 7 S with the inverse temperature f3 of the F 
model. Since the whole high-temperature phase is critical with central charge 
C — 1, in the thermodynamic limit 7 S should vanish for all (3 < (3 C — In 2, 
whereas in the non-critical ordered phase the exponent should stick to the pure 
quantum gravity value of 7 S = —1/2. To get an overview of the temperature 
dependence of 7 S we measured the distribution n^ 2 {B) of minBUs over an 
inverse temperature range of 0.2 < (3 < 1.4 for graphs of size N2 = 2048 
and performed fits of the functional form (19) to the data to extract 7 S [thus 
first neglecting the possibility of additional logarithmic corrections indicated 
in Eq. (20)]. The resulting estimates for 7 S presented in Fig. 8 show a plateau 
value of 7 S ~ —0.25 within the critical phase (3 < In 2 and a slow drop down to 
7 S pa —0.5 at (3 — 1.4 in the low-temperature phase. Note that the error bars 
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displayed in Fig. 8 are those resulting from the fit procedure itself and are 
thus not representing the full statistical variation due to the above mentioned 
cross-correlations between the values of n^ 2 {B). For comparison the correct 
error bars as obtained from a more elaborate jackknife analysis are shown for 
three selected /3-values, which are discussed in more detail below. As will be 
shown there, the fact that 7^ is found to be still considerably smaller than zero 
in the high-temperature phase is due to a finite-size effect. We do not employ 
the corrected fit (20) at this point, which is found to be unstable for the small 
graph size considered here. 

More precise estimates for 7,, are found from a FSS study of three series 
of simulations, one at the critical point (3 C = In 2, one in the critical high- 
temperature phase at (3 = 0.2 and one deep in the ordered phase at (3 = 1.4. 
For the latter case, the exponential slowing down of the combined link-flip and 
surgery dynamics of the graphs reported in Refs. j27, limited the maximum 



accessible graph size to N 2 = 16 384, while for the simulations at the critical 
point and in the high-temperature phase graphs with up to = 65 536 sites 
were considered. At (3 = 1.4, this maximal size is anyway sufficient since we 
find no finite-size drift in the estimate for 7,, with increasing graph sizes, all 
results being compatible with the conjectured value of 7 S = —1/2. Thus, as 
our final estimate for (3 = 1.4 we report the value found for N 2 = 16 384, 
7<- = —0.478(17). For the quoted statistical error estimates the jackknifing 
procedure described above for pure dynamical 4 graphs was used, thus taking 
full account of the present fluctuations. 

At the critical point (3 C = In 2 fits of the form (19) without logarithmic cor- 
rections show considerable finite-size effects, with 7 S slowly increasing with 
the graph size. For the largest graph size considered, N 2 = 65 536, the thus 
found estimate 7 S = —0.2075(17) is still far away from the expected result 
7 S = 0. Taking the logarithmic corrections into account, however, these re- 
sults can be considerably improved, with the numerical estimates for 7,, now 
being fully consistent with the theoretical prediction. The parameters of fits of 
the corresponding functional form (20) are collected in Table 5. Including this 
correction, no further finite-size dependence of the estimate 7 S is visible. The 
occurring values for the "correction exponent" k are not too far away from and 
indeed statistically compatible with the conjectured value of k = —2. Since for 
the case of N 2 = 65 536 only a much shorter time series than for the smaller 
graph sizes was recorded, we present as our final estimate of the critical value 
of 7, the result at N 2 = 32 768, 7, = 0.013(70). 

Finally, in the high-temperature phase at (3 = 0.2 the simulation results be- 
have very similar to the critical point case. When applying fits of the form 
(19) without logarithmic corrections, considerable finite-size effects are found, 
and the approach of the resulting exponent estimates 7^ to the expected value 
of 7 S = is very slow. On the other hand, the estimates resulting from fits of 
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Table 5 

Parameters of fits of the form (20) to the distribution n^ 2 {B) of minBUs for </> 4 
random graphs coupled to the F model at (3 = (3 C = In 2. 



N 2 






^7s 


7s 


K 




16 384 


100 


2048 


25.7(15) 


0.05(13) 


— 1.97(89) 


-10.9(69) 


32 768 


110 


4096 


27.08(93) 


0.013(70) 


— 1.80(50) 


-12.6(47) 


65 536 


120 


4096 


27.5(14) 


-0.05(12) 


-1.27(82) 


-6.9(71) 



the form (20) to the data are compatible with 7 S = for the larger of the con- 
sidered graph sizes. For N 2 = 32 768 we find 7, = -0.041(73), k = -1.38(47), 
Q = 0.05 with cut-offs B min = 100 and -B max = 2048. To complete the picture, 
it should be mentioned that the functional form (20) does not fit the data in 
the low-temperature phase at (3 = 1.4 well and does not give estimates of 7 S 
compatible with 7^ = 0, in agreement with theoretical expectations. 



4-3 The Haus dor ff dimension 



The internal Hausdorff dimension dh of the dynamical polygonifications is one 
of its most striking features. Apart from the physical implications, its large 
value causes a quite inconvenient obstacle for the numerical analysis of the 
model, namely the comparable smallness of the effective linear extent of the 
graphs at a given total volume N 2 as compared to flat lattices. As matter 
variables are coupled to the dynamical graphs, the strong coupling between 
graph and matter variables at criticality could lead to a change of the fractal 
dimension of the lattices. In a phenomenological scaling picture, such a strong 
coupling of matter and geometry should set in as soon as the correlation length 
of the matter system becomes comparable to the intrinsic length scale of the 
graphs or polygonifications. For conformal minimal matter, there has been 
quite some debate about how dh should depend on the central charge C of the 
coupled matter system, see, e.g., Refs. [26[ IH, H^, |66|, |6?J . For C = the 
result dh = 4 is exact [6^|. Furthermore, the branched polymer model |63| . 
describing the C — > 00 limit [l3|, yields dh = 2 (see, e.g., Ref. For the 

intermediate region < C < 1 two differing conjectures have been made for 
dh, namely 7(| 
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and UM 



d h = , ^ = ; ^ oo. (22) 



All numerical investigations up to now, on the other hand, are consistent with 
a constant d h = 4 for < C < 1 0, H Q- Naturally, the limiting 
case C = 1 considered here is of special interest for the investigation of the 
transition to the branched polymer regime C ^> 1. Numerically, it has proved 
exceptionally difficult to extract the Hausdorff dimensions from the statistics 
of the practically accessible graph sizes [26[ [73, EH ■ Only more recently, the 
development and application of suitable FSS techniques allowed for a more 
successful and precise determination of d h joll. I59I l65^. 



4-3.1 Scaling and the two-point function 

The fractal structure of the polygonifications is encoded in their geometrical 
two-point function. Here, different definitions are possible. While in Eq. (9) a 
definition in terms of the vertices of the graphs has been used, here, instead, 
the number of vertices of the quadrangulation is counted. Thus, we define the 
geometrical two-point function G^ 2 (r) as the average number of vertices of 
the polygonifications at a distance r from a marked vertex, where "distance" 
denotes the unique minimal number of links one has to traverse to connect both 
vertices. Since the intrinsic length of the model scales as N^ dh by definition of 
the internal Hausdorff dimension dh, from the usual FSS arguments one can 
make the following scaling ansatz (see, e.g., Ref. [65^). 

G$(r)~N?F(r/Ni /dh ), (23) 



i.e., Gjj(r) is a generalised homogeneous function and one can define a scal- 
ing function F(x) of the single scaling variable x = r/N^ dh and a critical 
exponent a. Due to the obvious constraint N 2 = E r (r), the exponent a 
is not independent, but given by a — 1 — It turns out that for practical 
purposes the scaling variable has to be shifted to yield reliable results, see, 
e.g., Refs. jH, 67, 73]. The necessity of such a shift can be most easily seen 



by a phenomenological scaling discussion of the mean extent defined by 

(r)iV 2 = ^E^n 2 W~^ 2 1/dh , (24) 
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with F = J2r F( r )- On general grounds, one expects the presence of analytical 
scaling corrections, 

i r )N 2 p _0_ _±_ (t) -, 
N l/dH ~ *° N l/d h N 2/ dh +■■■■ W 



Combining the terms proportional to 1/N^ dh on both sides, the mean extent is 
found to be (r + a)N 2 ~ F N 2 dh + 0(N 2 1 ^ dh ). Thus, to incorporate first-order 
corrections to scaling, the ansatz (23) is replaced by 

G^{r)~N«F[{r + a)/Nl /dh l (26) 
i.e., the scaling variable is now defined to be x = (r + a)/N^ dh . 
4-3.2 Scaling of the maxima 

The two-point function GfJ 2 (r) exhibits a peak at intermediate distances and 
declines exponentially as r — > oo, cf. Fig. 10 below. From the scaling ansatz 
(26) one infers the following leading scaling behaviour of the position and 
height of the maxima, 

''"max + a = A r N 2 h , (27) 
GS a W = A n Ni~ 1/dh + B n . 



Since the location and height of these maxima can be determined numerically 
from simulation data, these relations can be used to estimate the intrinsic 
Hausdorff dimension dh. A technical difficulty is given by the fact that r can 
only take on integer values for the discrete graphs considered. This problem 
is circumvented by a smoothing out of the vicinity of the maximum by a 
fit of a low-order polynomial to G^i(r) around its maximum. For practical 
purposes, we find a fourth-order polynomial sufficient for this fit. Reliable 
error estimates are found by jackknifing over this whole fitting procedure, 
where the individual statistical errors of the data points included in the fits 
are taken to be equal. Thus, one arrives at estimates for the peak locations r max 
and heights Gn(r max ) as a function of the graph size N 2 , to which then the 
functional forms of Eq. (27) are fitted. The effect of neglected FSS corrections 
is accounted for by successively dropping data points from the small- N 2 side. 
For simulations of pure 4 random graphs, in this way we find the value of dh to 
steadily increase on omitting more and more points. For the range N 2 = 4096 
up to N 2 = 32 768 we thus arrive at the estimates d h = 3.803(28), Q = 0.22 
from the scaling of the peak locations and d h = 3.814(63), Q = 0.44 from the 
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Fig. 9. Root mean square extent (r 2 ) 1 / 2 of regular c/> 4 random graphs with N% = 2048 
sites coupled to the F model. The horizontal line indicates the root mean square 
extent of pure </> 4 random graphs of the same size. 

peak heights. Both estimates are still noticeably away from the asymptotic 
values dh = 4, owing to the neglect of higher-order correction terms 
We note that introducing the shift parameter a already largely improved the 
estimates, since fixing a = we arrive at dh = 3.4313(20) from the peak 
locations. Further improvement is gained from the inclusion of the next-order 
correction term for the scaling of the peak locations, 



+ a 



A r N 2 1/dh 



+ B r N 2 ~ 



l/d h 



(28) 



which yields an estimate of dh = 3.964(42), Q = 0.24 for the range N 2 = 
512, . . . , 32 768, in perfect agreement with dh = 4. 

For random 4 graphs coupled to the F model, we find a small dependence 
of the root mean square extent on the inverse temperature (3 of the coupled 
F model and also a slight shift of (r 2 ) 1//2 as compared to the case of pure 4 
random graphs, cf. Fig. 9. Thus, one might expect the Hausdorff dimension dh 
to be temperature dependent, too. We performed simulations for three inverse 
temperatures, namely (3 = 0.2, (3 = In 2 and (3 = 1.4, covering the cases of 
interest. The results for dh from fits of the functional form (27) to the data are 
found to steadily increase on omitting more and more points from the small- 
N 2 side. In agreement with the case of pure 4 graphs, the final estimates 
for dh are found to be significantly smaller than dh = 4 for all three inverse 
temperatures due to the presence of higher-order corrections to scaling. For 
the peak heights, this analysis yields the estimates dh = 3.446(68) for (3 = 0.2, 
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Table 6 

Parameters of fits of the form (28) to the peak locations of the two-point functions of 
the random graph F model. The maximum graph size was N 2 = 65 536 for (3 = 0.2 
and = In 2 and N 2 = 32 768 for (3 = 1.4. 



p 




^4^> 


B r 


a 


dh 


Q 


0.2 


2048 


2.34(43) 


9.2(39) 


5.5(21) 


4.13(20) 


0.86 


In 2 


2048 


1.96(39) 


5.6(43) 


3.7(20) 


3.93(21) 


0.11 


1.4 


1024 


2.18(49) 


6.0(37) 


4.4(22) 


4.07(26) 


0.44 



d h = 3.426(92) for /3 = In 2 and d h = 3.94(23) for (3 = 1.4, where the rather 
different result for (3 — 1.4 again indicates the presence of competing local 
minima in the x 2 distribution. The found higher-order scaling corrections are 
resolved by using the fit ansatz (28) for the peak locations. Here, we do not find 
a significant sensitivity of the parameter estimates on the cut-off iV^min and 
for all three inverse temperatures the resulting values for dh are in agreement 
with the pure gravity value dh = 4, cf. the fit data collected in Table 6. 

4-3.3 Scaling of the mean extent 

As an alternative to the scaling of the maxima of the two-point function, one 
can also consider the behaviour of mean properties of the distribution G^ 2 (r), 
especially the scaling of the mean extent (24). Taking the next sub-leading 
analytic correction term into account, we make the scaling ansatz 

(r + a) N2 = A {r) N] /dh + B {r) N^ l/dh . (29) 

We again consider the case of pure 4 graphs first. When fixing Bi r \ = and 
adapting the lower cut-off A^min, the resulting values of dh are significantly too 
small in terms of the statistical errors with an obvious tendency to increase as 
more and more of the points from the small- N2 side are omitted. On the other 
hand, including the correction term of Eq. (29) largely reduces the dependency 
on the cut-off N 2Mn . For N 2Mn = 256 we find d h = 3.90(15), Q = 0.01, in 
nice agreement with d^ = 4. Here, the fits become very unstable as less points 
are included; this explains the use of the cut-off N 2)iain = 256, although the 
quality-of-fit is rather poor. 

The authors of Ref. [58] have proposed a different and less conventional method 
to extract a and dh from data of the mean extent, which they claim to be es- 
pecially well suited for obtaining high-precision results. They consider the 
combination R a ,N 2 (dh) = (r + a^iV^T 1 ^' 1 , and evaluate it for a series of sim- 
ulations for different graph sizes N 2 . Then, for a given a and for each pair 
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{N l 2 ,N 3 2 ) they define d%(a) such that R a ^ 2 {d L ^) = R aN j(d%), i.e., 
Ajf , lnJq-]nNi 

d-h ,{a) = . . : : ,, ■ r. (30) 

m((r)M + a)-m((r)^ + a) 



By a binning technique, an error estimate cr(d^) is evaluated and the estimates 
c?^(a) are averaged over all pairs (N 2 , N 2 ) of volumes, 4(a) = J2i<j d l h 3 (a)/Af, 
where Af denotes the number of pairs (N 2 , N 2 ). Then, the optimal choice a opt 
of the shift is found by minimising 

2/ \ \- K (a) -4(a)] 2 , . 

x(a) = § « ' (31) 



being accompanied by an optimal estimate 4(oopt)- The authors of Ref. |5£ 
suggest to estimate the statistical error of this final estimate by considering 
the variation of (a,dh) in an interval of a around a op t defined by x 2 ( a ) < 
min[l, 2x 2 (a op t)]- We implemented this procedure to compare with the results 
of the fits to Eq. (29) with = for the case of pure 4 graphs. We find 
the ad hoc assumption for the estimation of the errors of (a, dh) not adequate. 
Instead, we apply a second-order jackknifing technique (cf. Ref. ji^]) to be 
able to give error estimates for d%(a) as well as the final estimate (a, dh) which 
are found to be largely differing from those resulting from the rule x 2 { a ) < 
min[l, 2x 2 (a opt )], ranging from four times smaller to ten times larger error 
estimates. The estimates of dh itself are found to be indeed slightly increased as 
compared to the fit method (which yielded estimates clearly smaller than dh = 
4). This, however, can be traced back to the fact that the individual estimates 
d%(a) all receive the same weight in the average d~h(a) above, irrespective 
of their precision, giving an extra weight to the results for larger graphs, 
which cannot be justified on statistical grounds. If, instead, we use a variance- 
weighted average 



the resulting estimates for dh and a are statistically equivalent to those found 
from the fits to (29). For a cut-off A^ 2)min = 2048, for instance, we find dh = 
3.97(12) compared to dh = 3.99(12) from a simple fit of the form (29) with 
5( r ) = 0. Thus, we do not find any special benefits of this computationally 
rather demanding method as compared to a plain fit to (29) with = 
and hence do not present further detailed results for this method. 

For the case of the F model coupled to the 4 random graphs we proceeded 
as before, again using simulation data for (3 = 0.2, (3 — ln2 and (3 = 1.4. The 
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Table 7 

Parameters of fits of the form (29) including the correction term to the mean extent 
of dynamical (f) A graphs coupled to the F model at different inverse temperatures (3. 



p 


^2,min 


A {r) 


B {r) 


a 


dh 


Q 


0.2 


512 


2.58(48) 


11.4(33) 


7.0(22) 


4.08(21) 


0.10 


In 2 


512 


1.37(22) 


0.4(29) 


1.1(12) 


3.45(14) 


0.41 


1.4 


512 


2.6(10) 


9.1(58) 


6.2(42) 


4.15(47) 


0.29 



results from fits of the mean extent (r) n 2 to the form (29) with B/ r \ = show 
very much the same behaviour as the results from the scaling of the maxima 
of the two-point function, with estimates of dh clearly below dh = 4 and slowly 
increasing as more and more points from the small- N 2 side are omitted from 
the fits. The outcomes of the method of Ref. j^fj described above, with the 
average (32) and the x 2 { a ) TV ^ e replaced by a jackknife error estimate, are again 
very close to the fit results. Including the correction term of (29), i.e., relaxing 
the constraint Bi r ) = 0, on the other hand, yields estimates consistent with 
dh — 4 for (3 = 0.2 and (3 = 1.4, however with rather large statistical errors, cf. 
the parameters collected in Table 7. Note that, as mentioned before, the results 
for /3 = 1.4 are in general less precise than those for the other two inverse 
temperatures, which is due to the exponential slowing down of the combined 
link-flip and surgery dynamics in the low-temperature phase, cf. Ref. |27| . 
The fit for /3 = In 2 settles down at a completely different minimum of the 
X 2 distribution, yielding an almost unchanged dh compared to the outcome 
of the corresponding fit without correction term. This underlines the fact 
that the complexity of the chosen fit is at least at the verge of being too 
high for the available data. Nevertheless, combining the data for dh from the 
presented methods and including the comparison to the pure gravity case, we 
find no reason to assume that dh differs from dh = 4 for the case of the F 
model coupled to 4 random graphs. At any rate, the values dh ~ 4.83 and 
dh = oo resulting from the analytical conjectures Eqs. (21) and (22) for C = 1, 
respectively, are clearly incompatible with the results found here. 

Finally, we note that the parameters a and dh determined from the fits dis- 
cussed above lead to a nice scaling collapse of the two-point functions G^ 2 (r) 
when re-scaled according to the scaling ansatz of Eq. (26). Figure 10 shows 
this collapse of distributions for the case of (3 — 0.2 and the choice of pa- 
rameters found from a fit to the form (29) with B( r ) = 0, i.e., dh = 3.57(12) 
and a = 1.60(74). The visible deviations around the peaks of the distribu- 
tions indicate the presence of higher-order corrections not incorporated into 
the scaling ansatz (26). 
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Fig. 10. Scaling collapse of the two-point functions G^(r) of c/> 4 graphs coupled 
to the F model at /3 = 0.2, re-scaled according to Eq. (26) with dh = 3.57 and 
a = 1.60. 



5 Conclusions 



The six-vertex F model represents an example of the limiting case of a critical 
theory at the "barrier" of central charge C = 1, where the Liouville approach 
to Euclidean quantum gravity in two dimensions breaks down |12fl and the 
ensemble of planar random graphs coupled to such matter is at the verge of a 
presumable collapse towards a phase of minimally connected, tree-like surfaces 
termed "branched polymers" |l3j. At the same time, the family of ice-type 
vertex models of statistical mechanics includes as sub-classes a variety of well- 
known lattice spin models and combinatorial counting problems and hence the 
analysis of its coupling to two-dimensional Euclidean quantum gravity is that 
of a prototype model of statistical mechanics subject to annealed, correlated 
connectivity disorder from random graphs (see also Ref. jzfj). 

For studying the effect of a coupled six-vertex F model, we generalised the 
well-established methods of simulating dynamical triangulations to the case 
of planar quadrangulations and the dual "fat" 4 random graphs; the details 
of this simulational machinery will be presented in a forthcoming publication 
|27| . We have analysed the critical and off-critical behaviour of this model 
using a series of extensive Monte Carlo simulations and subsequent finite-size 
and thermal scaling analyses. On the square lattice, this model undergoes 
an infinite-order phase transition of the Kosterlitz-Thouless type to an anti- 
ferroelectric phase of staggered order. Expecting similar ordering behaviour 
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to occur for the model on a random quadrangulation, we generalised the cor- 
responding staggered polarisation (the order parameter) to the random graph 
case by a duality transformation of the vertex model. 

The scaling analysis of the simulation data is hampered by the presence of 
extraordinarily strong corrections, which can be traced back to the combined 
effect of the comparable smallness of the effective linear extents of the con- 
sidered (two-dimensional) lattices due to their large fractal dimension close 
to four and the presence of logarithmic corrections generically expected for 
a C = 1 critical point. Additionally, the form of the critical singularities for 
a Kosterlitz-Thouless phase transition severely limits the effectivity of the 
usual finite-size scaling techniques. General symmetry considerations imply 
that the C = 1 critical point of the random-graph model should occur at the 
coupling (3 C = In 2, which is quite remarkably identical to the critical coupling 
of the square-lattice model. Additionally, this is in agreement with a matrix 



model treatment of the system [17|,|20|. Due to the aforementioned strength 
of scaling corrections, a precise determination of the critical coupling from the 
scaling of the polarisability alone is found to be hard. A comparison of the 
peak positions re-scaled according to the mean linear extents of the lattices 
between the random graph and square-lattice models 0|, however, shows that 
the finite-size scaling approaches of both models are indeed very similar, but 
with larger correction amplitudes for the random graph model. Subsequently, 
however, a precise and consistent estimate of the transition point could be 
extracted from the scaling of the co-ordination number distribution of the 
graphs. A cursory comparison of the scaling behaviour of the model for differ- 
ent ensembles regarding the inclusion of singular contributions in the graphs 
reveals that corrections to scaling increase as more and more singular con- 
tributions are included. This contrasts with the findin gs fo r pure gravity and 
Potts models coupled to the polygonifications model |49j, |56( . As far as the 
critical exponents related to the order parameter are concerned, a finite-size 
scaling analysis of the spontaneous polarisation and the polarisability at the 
asymptotic critical coupling yields critical exponents in good agreement with 
the predictions from the KPZ/DDK formula. An attempted thermal scaling 
analysis of the polarisability around its peak remains inconclusive due to the 
extraordinary magnitude of finite-size corrections. As a curiosity, we report 
the finding of a critical internal energy of the model, U(/3 C ) = 1/3, which is 
identical between the square-lattice and random graph cases. 

Several aspects of the back-reaction of the matter variables onto the properties 
of the 4 random graphs are analysed as a function of temperature. The dis- 
tribution of co-ordination numbers of the quadr angulations can be determined 
very accurately. The fraction of quadrangulation sites of co-ordination num- 
ber two is found to be peaked around the asymptotic critical coupling, thus 
defining a pseudo-critical point which determines the infinite-volume critical 
coupling very accurately. A scaling analysis of the distribution of "baby uni- 
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verses" of the graphs in the spirit of Refs. [61, 64] allows to extract the string 
susceptibility exponent 7 S of the model. It is found to coincide with the value 
7 S = expected for a C = 1 theory throughout the critical high-temperature 
phase. The pure-gravity value 7 S = —1/2 is found in the non-critical low- 
temperature phase. Exploiting finite-size scaling relations, we finally analyse 
the geometrical two-point function of the graphs and extract the fractal Haus- 
dorff dimension. We find it to be consistent with the pure gravity value dh — 4 
for all temperatures of the coupled vertex model. The analogous analyses for 
the case of pure 4 random graphs convincingly demonstrate the universal- 
ity of these graph-related critical exponents with respect to a change from 
triangulations to quadrangulations. 

In summary, despite of the presence of scaling corrections of extraordinary 
size, a careful analysis of our simulation data allows for an independent con- 
firmation of the location of the critical point and the behaviour of the string 



susceptibility exponent predicted by the matrix model treatment [17|, |20j . In 
addition, the behaviour and critical exponents related to the order parame- 
ter, energy-related observables, as well as further geometrical properties such 
as the Hausdorff dimension can be reliably determined. An even richer be- 
haviour can be expected for the 8-vertex model coupled to dynamical quad- 
rangulations, such that its analysis by a series of simulations similar to the 
one presented here would be a promising future enterprise. 
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